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o3 " Abstract 

i— ), 

t^. ■ A new analytic treatment of the two-dimensional Hubbard model at finite temperature 

and chemical potential is presented. A next nearest neighbor hopping term of strength t' is 

QJ \ included. This analysis is based upon a formulation of the statistical mechanics of particles 

"^ | in terms of the S-matrix. We show that for U/t large enough, a region of attractive 

interactions exists near the Fermi surface due to multi-loop quantum corrections. For 

g; 

i ■ t' = —0.3, these attractive interactions exist for U/t > 6.4. Our analysis suggests that 

superconductivity may not exist for t = 0. Based on the existence of solutions of the 

o' 

integral equation for the pseudo-energy, we provide evidence for a phase transition and 

(SI ■ estimate T c /t « 0.02 for U/t = 7.5 and t'/t = -0.3 at hole doping 0.15. 
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I. INTRODUCTION 



The Hubbard model in two spatial dimensions plays a central role in the modern 
theory of strongly correlated electrons. Since it is believed to be a good microscopic 
model for the underlying physics of high T c superconductivity (HTSC) jl|, it has been 
studied extensively over the past two decades. A partial list of publications on the 



thermodynamics of the Hubbard model is [2|-|l0l|. Despite this effort, many of its 
important properties are not currently well-understood, and it remains to be estab- 
lished definitively that it possesses all of the main features of HTSC For instance, the 
precise mechanism that leads to attractive d-wave pairing, which must arise from the 
purely electronic Coulomb repulsion, is still not well understood. Analytic methods 
to date are rather limited since the model is in the strong coupling regime. Lattice 
Monte-Carlo methods on the other hand are limited to small lattices and suffer from 
the fermion sign problem, especially at non-zero doping and sufficiently low temper- 
atures. For these reasons, any new analytic methods, though approximate, may shed 
new light on the problem. 

In this work we present an analytic approach to the thermodynamics of particles at 
finite density and temperature based on the reformulation of the statistical mechanics 
of particles in terms of the zero temperature and density particle-particle S-matrix 



developed in [111 . Il2j . As explained there, the potential advantage of this method is 
that, unlike the usual diagramatic Matsubara approach to finite temperature field 
theory, it disentangles the zero temperature dynamics from the quantum statistical 
sums. Some remarks clarifying the nature of this formalism, and the approximations 
made, are called for. The approach was modeled after the thermodynamic Bethe 
ansatz 13|, which is exact for integrable models in one spatial dimension since the N- 
body S-matrix factorizes into 2-body S-matrices, Our method indeed reduces to the 
thermodynamic Bethe ansatz in the two-body scattering approximation, as shown 



in [12| for the interacting Id Bose gas. The main approximation we make is that 
we consider only many-body processes that involve arbitrary numbers of primitive 
binary collisions. I.e. we neglect processes that in some sense involve 3 or more 
particles colliding simultaneously, which is not the same as ignoring the many-body 
aspect of the problem altogether. In this non-relativistic context, it is well-known 
that the two-body S-matrix can be calculated exactly, thus in some regards the 
method is non-perturbative. Although this is not a fully controlled approximation, 
it has been demonstrated to give reliable results for other strong-coupling problems, 
in particular the critical point of the 2d bosonic gas and more importantly to Bose 
and Fermi gases in the 3d scale- invariant unitary limit |l2J Il4j . For instance for 
the unitary Fermi gas on the BEC side of the crossover, the critical temperature 
calculated is consistent with Monte-Carlo methods, and the ratio of the viscosity 
to entropy density agrees very well with the most recent experimental data[15|. It 
should be pointed out that the exact Bethe-ansatz solution of the Id Hubbard model 



exhibits two additional holon excitations 16[ , which can be inferred from poles in the 



S-matrix of the fundamental fermions; we have no evidence for such excitations in 
the 2-dimensional case, thus it is unclear whether a meaningful comparison with the 
Id case can be made. 

Our conventions for the Hubbard model are described in the next section. We 
include a next nearest neighbor hopping term of strength t', since it is known to 
be non-zero in the cuprates; as we will show, its effects are important. In section 
III, the effective momentum dependent coupling, which is the kernel G(ki,k2) con- 
structed from the logarithm of the 2-body S-matrix in the integral equation satisfied 
by the pseudo-energy, is analyzed. We show that there exists a band of attractive 
interactions near the half-filled Fermi surface for U/t large enough, and t' plays an 
important role in determining this property. We emphasize that no approximations 
are made in section III, since, as stated above, the 2-body S-matrix can be calculated 



exactly, and this attractive region exists regardless of the subsequent approximations 
we make in studying the thermodynamics. This attractive mechanism, which arises 
from quantum loop corrections, appears to be different than other mechanisms dis- 
cussed in this context, such as those based on spin fluctuations or the resonating 
valence bond picture. It is thus important to investigate the consequences of these 
attractive interactions and how they might be connected to HTSC, and this paper is 
a first step in this direction. In section IV the S-matrix based formalism we utilize 
for calculating thermodynamic properties is reviewed and specialized to the Hubbard 
gas. In section V the free energy is analyzed, and we present some evidence for phase 
transitions. 

II. HUBBARD MODEL CONVENTIONS 

The Hubbard model describes fermionic particles with spin, hopping between the 
sites of a square lattice, subject to strong local coulombic repulsion. The lattice 
hamiltonian is 

H=~t ^ ( C k a Cr j>a ) ~t' Yl { C L a C r 3 ,a) +UY n rt U rl (1) 

<ij>,a=t4 <ij>',a=t,4- r 

where rjj,r are sites of the lattice, < i,j > denotes nearest neighbors, n = c^c are 
densities, and c^,c satisfy canonical anti-commutation relations. For both cuprates 
LSCO and BSCO, U/t ~ 13. We have also included a next to nearest neighbor 
hopping term t', since it is not difficult to incorporate into the formalism, and it is 
known to be non-zero for high T c materials. As we will see, it can play a significant 
role. For LSCO and BSCO, t'/t approximately equals —0.1 and —0.3 respectively; 
in our analysis below we set t'/t = —0.3. 



We introduce the two fields ipt,l an d the action 

(2) 



s= fdhdtl J^^i^a-n) 

J Va=u / 



where % is the hamiltonian density. The field has the following expansion charac- 
teristic of a non-relativistic theory since it only involves annihilation operators, 



d 2 k 



Mr)= / ^c k , a e Jk - r (3) 

J 2tx 

and satisfies 

{Mr)^Ur')} = 5(r-r')5 ata/ (4) 

Since we have represented sums over lattice sites r as J d 2 r/a 2 , where a is the lattice 
spacing, c r = aip(r). The free part of the hamiltonian is then 

H iicc = / d 2 k w k ^2 c l,a c K» ( 5 ) 

with the 1-particle energy 

Wk = — 2t (cos(k x a) + cos^a)) — At' cos(k x a) cos(k y a) (6) 

where t taken to be positive. In the sequel it is implicit that k is restricted to the 
first Brillouin zone, —tt/cl < k XtV < n /a 

The interaction part of the hamiltonian is local, and becomes a continuum integral: 

H int = ^Jdh VfrM>M (7) 

where w = 2Ua 2 . The model is now viewed as a quantum fermionic gas, where the 
only effect of the lattice is in the free particle energies u^. 

The field ip has dimensions of inverse length, and the coupling u units of energy • 
length . In the sequel we will scale out the dependence on t and the lattice spacing 
a, and physical quantities will then depend on the dimensionless coupling 

u 2U 

Positive g corresponds to repulsive interactions. 



III. THE EFFECTIVE MOMENTUM DEPENDENT COUPLING AND 
THE POSSIBLE ORIGIN OF ATTRACTIVE INTERACTIONS 



In the finite temperature formalism developed in ll|, |l2] , the occupation numbers 
/ are parameterized in terms of a pseudo-energy e(k) in the same manner as for a free 
theory: / = \j[eP e + 1). In the approximation that only the many-body processes 
built out of primitive 2-body collisions are retained, the pseudo-energy satisfies an 
integral equation based on a kernel G(ki,h 2 ) which is related to the logarithm of 
the 2-particle S-matrix. This approach to the thermodynamics will be reviewed in 



the next section. The final result derived in [12[ involves only the particle-particle 
S-matrix at zero temperature and density, which can be calculated exactly. The 
temperature and density enter the formalism in the integral equation for the pseudo- 
energy, thus this formalism does not require particle-particle or particle-hole Green's 
functions at finite temperature and chemical potential. In this section we study the 
main features of the kernel and demonstrate that there are regions of the Brillouin 
zone where the interactions are effectively attractive, even though the bare model 
has repulsive interactions. We we wish to emphasize that no approximations are 
made in obtaining the results presented in this section, which essentially amount 
to quantum corrections to scattering, and some conclusions are independent of the 
thermodynamics studied in subsequent sections. 

A. Structure of the kernel 



As described in |l2J, the kernel has the following structure: 



G(k 1 ,k 2 ) = -^log(l + 2zZM) (9) 



where Ai is the 2-body scattering amplitude. (We are suppressing the momentum 
dependence.) X represents the available phase space for two-body scattering: 

1 = h I d2p 5{ ~ E ~ Up ~ Wk " p) (10) 

where E and K are the total energy and momentum of the two incoming particles 
with momenta ki and k^: 

E = co kl +uj k2 , K = ki + k 2 (11) 

All energy and temperature scales, E, T, the chemical potential /i, and t', will 
be expressed in units of the hopping parameter t. We thus scale t out of u^ so 
that henceforth Wk equals (EI) divided by t. We will also rescale k by l/o so that 
— 7r < k xy < 71". Since X is proportional to 1/t, G oc t, and henceforth G will represent 
G/t, which is dimensionless after scaling out factors of a also. The kernel G then 
depends only on the dimensionless coupling g defined in eq. ([8]), and the momenta. 

The scattering amplitude can be computed by summing multi-loop ladder 
diagrams [12|], which factorize into 1-loop integrals in this non-relativistic context. 



This leads to 



M = - 9/ \ (12) 

1 + igL/2 v ; 



where g is the coupling defined in eq. ([8]). L is a 1-loop integral: 

dcu f d 2 p 



2tt J (2ir) 2 \uj — lo p + iej \E — u — w K - P + ^ e 
d 2 p 1 



(13) 



(27r) 2 E — Up — uk-p + 2ie 

where e is small and positive. In the numerical analysis below we set e = 0.001. 
Using Im(x + ie) _1 = — ir5(x), one sees that 

L=X + ij (14) 



where X is the phase space factor in eq. (fTUj) and is real and positive, and 7 is its 
imaginary part. Putting all of this together one has 

The imaginary part of the loop integral renormalizes the coupling g 

9S = T^f2 < 16 > 

Note that the manner in which gn enters the kernel leads to a well-defined large 
coupling limit; this was exploited for unitary quantum gases in [14J where gn is 
proportional to the scattering length which goes to ±00 in the unitary limit. 

The argument of the log in eq. ( 1151) can be identified as the 2-body S-matrix, which 
is unitary, i.e. S*S — 1. It should be emphasized that this is the exact two-body 
S-matrix, and this is possible because the model is non-relativistic. More specifically, 
the fields ip^^ only involve annihilation operators, in contrast to relativistic theories 
which are expanded in both creation and annihilation operators. Let us elaborate on 
this important point, which is well-known in other contexts, such as non-relativistic 
quantum gases, but often not completely explained. Represent the interaction vertex 
with two incoming arrows for the annihilation operator fields ip^^ and two outgoing 
arrows for the creation fields ipli- Consider for simplicity the 1-loop contributions. 
Diagrams with a closed loop, i.e. with arrows circulating in the same directions, 
such as the second diagram in Figure dj are zero because the integration over energy 
u inside the loop has poles in the integrand that are either both in the upper or 
lower half-plane, so that the contour can be closed at infinity without picking up 
residues. (This is only true because our formalism only involves the particle-particle 
S-matrix at zero temperature and density.) In other words, there is no "crossing- 
symmetry" as in relativistic theories, where there are three non-zero 1-loop diagrams, 
with different momentum dependence, which are crossed versions of each other. The 

8 



non-zero multi-loop diagrams are only of the "ladder type" , which factorize, and this 
was implicitly used in 12|. There is actually no fermionic minus sign associated with 
each loop since the arrows do not form a closed loop. Since this S-matrix is exact 
to all orders in g, the thermodynamic formalism we will use embodies some non- 
perturbative aspects of the problem, although it still represents an approximation to 
the thermodynamics, as explained in more detail in the next section. 



=0 





FIG. 1: One-loop contributions to the S-matrix. Only the diagram to the left is non-zero. 



B. Origin of attractive interactions 



The kernel G by construction is real. For small coupling g, G is independent of 
momentum and equal to —g/2. Thus G may be viewed as an effective, momentum- 
dependent coupling constant, and provides valuable information on the effective 2- 
body interactions at zero temperature. When G is negative the interactions are 
effectively repulsive, otherwise they are attractive. The important point is that the 
renormalization of g to g^ can in fact change the sign of g^, which changes the sign 
of G due to the branch cut in the logarithm. To demonstrate how this can happen, 
we first perform the p y integral in L. The result is 



2tt 2 



dp x 



1 



VB 2 + C 2 -D 2 



log 



(C-D) 



VB 2 + C 2 -D- 



log 



(D-C) 



VB 2 + C 2 -D 2 

(17) 



where 

B = 2 sin K y + At' cos(K x — p x ) sin K y 

C = 2 + 2 cos A' y + At' cos p x + At' cos(K x - p x ) cos K y (18) 

D = E + 2 cosp x + 2 cos(-R' a ; — p x ) + lit 

(Recall t' represents t'/t.) The sum of the logarithms in the above formula is simply 
log(— 1) = ±27r; however expressing the integral in this fashion ensures one is on 
the proper branch. The above formula proved to be very useful for the numerical 
evaluation of the kernel. 

From the expression for the renormalized coupling g^, eq. ( TT6j) . one sees that g^ 
can become negative if 7 is positive and g large enough, g > 2/7. Let ki, hi be in the 
center of mass frame, \i\ = — k 2 = k, so that K = 0, and L only depends on the total 
energy E. The loop integral 017|) can be expressed in terms of elliptic functions. To 
regulate the integral, we let the upper limit of the p x integral be it — n, and then let 
k — > 0. The loop integral L oc F(i log(4a//«), b) where F is the elliptic integral of the 



first kind, with a = y/(E - 8 + 8t')/(E - W) and b = (E - 8t') 2 /((E + 8t') 2 - 64). 
As k — > 0, we use lim x _^ 00 F(ix, b) = iK(l — b), where K is the complete elliptic 
integral of the first kind. The final result is 

r 2/ E-8t> X' 2 K ( 32(^-2) \ 

K "° n\(E + 8 + 8t')(E(8-E) + 6At'(t'-l))J \(E + 8f) 2 - 6Aj { ' 

with E — > E + lie. The flip in sign is a result of the combination of logarithms in 
eq. (HZ). 

The imaginary part of L, i.e. 7, is plotted in Figure [2] for t' = —0.3. Since 
7 is positive for large enough E, one reaches the remarkable conclusion that for g 
large enough, the effective interactions can become attractive. One can estimate 
this threshold for g as 2/7 max , where 7 max is the maximum value of 7 which occurs 
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where 7 flips sign. This value of 7 max can be obtained using the formula (TT9l) . and 
is 7max ~ 0.156 for t' = —0.3. This minimal threshold in g should be contrasted 
with the Cooper instability, which leads to superconductivity for arbitarily weak 
coupling. This fact is a consequence of a logarithmic divergence in the analogous 
loop integrals, which leads a gap proportional to e~ x l u for some appropriate coupling 
u. In the present context, the attractive interactions arising from the above change 
of sign do not involve an analogous logarithmic divergence, hence the threshold. 

In summary, we have shown that there are effectively attractive interactions above 
a threshold in U/t, e.g. for t' = —0.3, attractive interactions exist for g > 12.8, i.e. 
U/t > 6.4. We repeated this analysis for other values of t' using the formula (I19J1 . 
and our results for the mininum value of g necessary for attractive interactions are 
shown in the table below. A minimal threshold for superconductivity was proposed 
in (l7|, with minimal values of U/t in the comparable range of 4 — 7. On the other 
hand, the study in 18] indicates no threshold, namely, for the particular mechanism 
they they study, superconductivity exists for arbitrarily low U/t. This does not 
necessarily contradict our result, since the attractive mechanism we study here is 
essentially different, and we have not made a case yet that it is the one responsible 
for superconductivity. 



t'/t 


fi'min 


-0.1 
-0.2 
-0.3 
-0.4 


16.7 
14.3 
12.8 
10.0 



The change in sign of the effective coupling described above is reminiscent of 
what is encountered in the BEC/BCS crossover of the non-relativist ic continuum 
three-dimensional unitary gas. The model is defined with repulsive interactions, 



11 



-0.2 



-0.4 



-0.6 



-0.8 



E/t 



FIG. 2: The imaginary part of the loop integral 7 as a function of energy E for t'/t = —0.3. 



i.e. positive coupling, however there is a fixed point at a negative coupling g*, 
independent of momenta. The scattering length is proportional to the renormalized 
coupling gn = g/{l — g/g*). Just above the fixed point, the scattering length goes to 
—00, whereas just below it goes to +00. The kernel G for this model flips sign as one 
crosses the fixed point for the same reasons as above, i.e. because of the branch-cut 
of the logarithm 14j. Thus, the effective interactions can be repulsive or attractive, 
depending on which side of the fixed point one sits, even though the bare model 
defined by the hamiltonian had only repulsive interactions. The main difference in 
the Hubbard gas is that the renormalized coupling depends on the momenta, so that 
the interactions may become attractive in distinct regions of the Brillouin zone. 

We mention that a change in sign of certain couplings under renormalization group 
flow was found for a 2-chain Hubbard model (2-legged ladder) in 19j. Such ladders 
are effectively 1 dimensional, and were mapped onto an anisotropic Gross-Neveu 
model, i.e. free Dirac fermions with marginal current-current interactions. Since 
these Gross-Neveu models are very different from those considered here, it seems 
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unlikely that the change of sign described in 19| is related to the one described here, 
since our model is intrinsically 2-dimensional, and the phenomenon does not involve 
any renormalization group flow. 

We now study the kernel G(ki, k^) for ki = — k2 = k and verify the above results. 
For these momenta, the kernel is only a function of the total energy E. In Figure 
[H we plot G(E) for the values of the coupling g = 5, 13.5, 14, 15, 20 and if = —0.3. 
One observes that for the smaller g = 5, the kernel is everywhere negative. One can 
verify that for g large enough, the main features do not depend strongly on g. The 
most interesting region of g is around g ^ 13 — 15 for t' = —0.3. Comparing g = 14 
and 15, one sees that for g = 14 the attractive band is narrower, and for g — 13 
essentially disappears. We will fix g = 15 in our subsequent thermodynamic analysis, 
since this is in the interesting region and the attractive band is not too narrow. 

G 




FIG. 3: The kernel G as a function of total energy E for g = 5, 13.5, 14, 15, 20, and 
t' /t = —0.3 (Color figures on-line.) 

Figure H] shows the kernel for g — 15 for if /t = 0, —0.1, —0.3, —0.4. One sees 
that the interactions are effectively attractive only around a small region centered 
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at E R3 —1.5 for t'/t = —0.3. An important feature of our analysis is that it clearly 
shows the importance of a non-zero t'. For fixed g, if \t'\ is too small, there is 
no attractive region, as is apparent in Figure HI If superconductivity indeed arises 
from these attractive interactions, then this suggests that superconductivity may not 
exist if t' = 0. There is actually some experimental evidence for this, in that T c as a 
function of t'/t appears to extrapolate to zero[20| . It should be pointed out however 
that as one lowers t'/t, attractive regions continue to exist as long as one raises the 
coupling g, as is evident in the table above for g min as a function of t'. 

G 
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FIG. 4: The kernel G as a function of total energy E for g = 15 and t' = 0, —0.1, —0.3, —0.4. 

It is also intructive to plot G for pairs of opposite momentum as a function of 
k x , k y in the first Brillouin zone. This is shown in Figure [5] Again this shows that 
the interactions are attractive in a narrow region around half-filling. The positive 
regions at the corners of the Brillouin zone are due to a divergence in the loop integral 
which should be regularized; however since we will be studying hole doping of the 
half-filled state, the densities will be low enough to be far from these regions, so 
this regularization will be unnecessary. Finally, note that for low enough E, the 
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interactions are always repulsive, which should imply that at high enough doping 
the theory should be well-approximated by a Fermi liquid. 




FIG. 5: The kernel G for Cooper pairs in the first Brillouin zone for g = 15 and t'/t = —0.3. 
The horizontal axes are — n < k x , y < n and the vertical axis is the effective coupling G. 

If the attractive interactions exist near the Fermi surface, then Cooper's original 
argument should apply: the filled Fermi sea just serves to block states and the 
particles can form a bound state, i.e. Cooper pairs 2l|. Let us then make the 
hypothesis that the regions of attractive interactions described above lead to Cooper 
pairing. Then the following scenario emerges. The Fermi surface with interactions 
is calculated in section V based on the filling fractions /. In Figure [6] we plot these 
Fermi surfaces for various hole doping, and also display the attractive band. (See the 
next section for the precise definition of hole doping h; as defined it corresponds to 
the number of holes per plaquette.) These computed Fermi surfaces closely parallel 
experimental measurements, in that they flare out in the anti-nodal directions i .e. 
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(k x ,k y ) = (0,7r) and 90° rotations thereof 22l| . This figure shows that the attractive 



regions in the anti-nodal directions play the most significant role. At low densities 
(high hole-doping), there are no attractive interactions within the Fermi surface, 
and the model should correspond to a Fermi-liquid. On the other hand, as h is 
decreased, the Fermi surface intersects attractive regions in the anti-nodal directions. 
This first occurs around a hole doping h = 0.3. If a gap forms in these directions, 
then this could explain the anisotropy of the gap, which is zero in the nodal (it, it) 
directions. As the density is increased further, eventually the Fermi surface is beyond 
the attractive band. i.e. the attractive band is completely enclosed by the Fermi 
surface. At half filing, h = 0, the attractive regions in the anti-nodal directions are 
just inside the Fermi surface, however the Fermi surface still intersects the attractive 
band in the nodal, (it, it), directions. 

IV. THERMODYNAMICS FROM THE S-MATRIX 

The important question that remains is whether our approach to the thermody- 
namics of the Hubbard model can capture the instabilities proposed at the end of 
the last section based on the properties of the kernel G. If so, this will provide a 
calculation of the critical temperature. 

In this section we describe how to compute the free energy from a formalism that 



is a synthesis of the works [111. Il2j . Being a synthesis, it is worthwhile reviewing 
the main features of the derivation, and how the construction follows from the basic 
ingredients in these two papers. 

The starting point is a formal expression for the partition function Z in terms of 
the S-matrix derived in 



z ° + h I dEe ~ PE lmdE log ^ (20) 



16 




FIG. 6: Fermi surfaces for various hole doping h = 0,0.1,0.2,0.3,-04, as computed in 
section V. The pink region (grey offline) is the band of attractive interactions for g = 15. 
The axes are k x and k y in the first Brillouin zone, i.e. in the range — it to ir 

where S(E) is the off-shell S-matrix operator, Zq the free partition function, and 
/3 = 1/T. Although the above formula is simple enough, a considerable amount of 
additional work is needed to obtain something useful out of it. For instance the clus- 
ter decomposition property of the S-matrix is needed to show that Z exponentiates 
to an extensive free energy. 
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Consider for simplicity a single species of fermions. The basic dynamical variables 
are the occupation numbers / which determine the density: 

n -lw? m (21) 

Using a Legendre transformation in the variables n and /i, where /i is the chemical 
potential, one can show that there exists a functional F(f) such that the physical 
free energy follows from the variational principle SF/Sf = 0. This functional can be 
separated into a free part Fo and an interacting part Fi, 

F = Fo + Fi (22) 

The interacting part contains contributions from N to N particle scattering for all 
N. One expects the 2-particle term to be the most important and is of the form: 

F ^-\Sw?Sw? mG()iM)m (23) 

Although G is built only from the 2-body S-matrix, the formalism re-sums all many- 



body processes that involve arbitrary numbers of primitive binary collisions. (In ll| . 
certain terms in fTSOj) (referred to as Zb terms) were incorrectly dropped. This was 
corrected in [12j, which led to the expression in the last section for the kernel G.) 



The primary difference of the two works ll|, [l2j is the choice of Fo • The choice 
made in [12| was better suited to the diagrammatic expansion, and the resulting 
integral equation effectively sums up an infinite number of diagrams. However it 
was found for the present problem that this integral equation only has solutions in 
a very limited range of temperature and chemical potential, indicating that the sum 
of diagrams does not converge. In contrast, it turns out the choice of Fo made in 



111 ] does not suffer from this problem. The latter Fo also has an appealing physical 



interpretation, as we now explain. Consider 



ro = /^((c*-A*)/-^[(/-l)log(l-/)-/tog/]) (24) 
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where u^ is the 1-particle energy of the free theory. The above expression can be 
interpreted as Fo = e— Ts, where e is the first (u;— //)/ term and represents the energy 
density The remaining term represents the entropy density s|24j. This choice of Fo 



also more closely parallels the derivation of the thermodynamic Bethe ansatz 13) . 
Let us parameterize the occupation numbers in terms of a pseudo-energy e: 

m = isthr (25) 

Then the variational equation SF/Sf = can be expressed in the simpler form: 

e(k) =tu*-»-tJ-^ G(k, k') ^ - i (26) 

(We have restored the hopping coupling t here.) Using the above equation in F, the 
free energy density T can be expressed as 



d 2 k 
7 =-T 



(2 



71 



log(l + e~ Pe ) + y^-^z - ^k + /i) 



(27) 



Comparing with [12|, one sees that in the limit of small G, the equations presented 
there reduce to eqns. (j26)27p . 

For two-component fermions, the occupation numbers are parameterized in terms 
of two pseudo-energies e-y,, and they satisfy a coupled system of two integral equa- 
tions: 

e t (k) =tUk -^-tJ^ G(k, k') -^1 - i (28) 

and the same equation with t^^i- Here the kernel G is related to the scattering 
of spin up with spin down particles. By the SU(2) symmetry, for equal chemical 
potentials /i-j- = jjl^ = /i, 5^ = £^ = e, and one only needs to solve one integral 
equation. The occupation number for each spin component has the form of a free 
theory given in ( 1231) . and the total density is 2 times the expression in ( 12 ip . as is the 
free energy. 
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For the Hubbard model, since we scaled out t and the lattice spacing, everything 
then depends on the dimensionless variables jw = fi/t, T = T/t, v = t'/t. Note that 
all the temperature dependence is in T, thus possible phase transitions should occur 
at fixed values of T/t, for given g, t'. Henceforth we drop the hats, it being implicit 
that T, fi and t' are in units of t. The free energy density then takes the form: 

F=-L c (fi,T) (29) 



where we have defined a scaling function c (we suppressed the dependence on t')\ 
f d 2 k 

It will also be convenient to express the density n = —dT/dfi as 



log (1 + e"^)) + I ^T^(k) - u; k + M ; 



(30) 



n= M!g)_ = lz£ (3!) 

a 2 a 1 

where 

q ~ J J2vf e^) + 1 (32) 

Since 2q is the number of particles of either spin per lattice site, half-filling cor- 
responds to q = 1/2. The quantity h then corresponds to hole doping when it is 
positive, otherwise it represents particle doping. More precisely, h is the number of 
holes per plaquette and the lattice is completely depopulated at h = 1. 

In order to probe the properties of the model, we will need a few other thermo- 
dynamic quantities. As usual the pressure p = —J 7 . Consider first the entropy per 
particle, S/N = s/n, where the entropy density s = —dTjdT. It can be expressed 
in terms of the scaling functions as follows: 

§=l q{ c + Tfrc) (33) 

The energy density e = E/V = Ts + fin + J 7 . Thus the energy per particle is 

F T 2 

— = — = fi+—c (34) 

Nt nt P 2q K J 
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The specific heat per particle Cy/N at constant volume and particle number N 
is slightly more complicated since one must impose the constant density constraint. 
Setting dq/dT = relates \x and T derivatives as follows: 



dfi drq 

Using this, the specific heat per particle has the following expression: 



(35) 



C v 1 (dE\ T ( T (drqT 



V. THE FREE ENERGY AND ESTIMATES OF CRITICAL TEMPERA- 
TURES 

In this section we analyze the thermodynamics based on the formulas of the last 
section, provide evidence for instabilities, which may perhaps be phase transitions, 
and estimate critical temperatures. 

In order to study the free energy, one must first solve the integral equation ( J26l) 
for the pseudo-energy e. This can be done iteratively, i.e. one starts with the 
approximation Eq = u^ — ft and plugs this into the right hand side to generate e\\ 
this procedure is repeated until the solution converges. We approximated the integral 
equation by approximating the Brillouin zone as a 10 x 10 grid, and performing the 
integrals as discrete sums. This is rather crude, and was due to our limited computing 
resources; certainly one can do better. It was found that for large portions of the 
/i, T parameter space, the iterative procedure converged rapidly, typically within 10 
iterations. 

For reasons stated above, our analysis was performed for g = 15, and t'/t = —0.3. 
For fixed doping h, the chemical potential depends on temperature, but we find 
this dependence to be weak. In Figure [7] we plot hole doping ft as a function of 
chemical potential at the low temperature T = 0.2, and it is nearly a straight line. 
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As expected, increased doping corresponds to decreasing chemical potential; half 
filling occurs around fi = 0.5. This positive value of \i at half-filling is due to the 
mainly repulsive interactions. 

h 




-2.0 -1.5 -1.0 -0.5 



/' 



FIG. 7: Hole doping h as a function of chemical potential for T = 0.2. 

The computed Fermi surfaces for various hole doping are shown in Figure El In 
this figure they are defined as the contour where the filling fraction / = 1/2 at 
the low temperature T = 0.025. They are in good agreement with experiments 221 ]. 
especially in the anti-nodal directions, whereas in the nodal directions they are pulled 
back toward the center in a more pronounced manner in the data. This can likely 
be accounted for by adding additional hopping terms, such as next-to-next nearest 
neighbor. 

The most interesting feature of the integral equation (f26l) is that there are regions 
in /i, T where the iterative procedure does not lead to a solution for an arbitarily 
high number of iterations; the procedure leads to e that jumps successively between 
two values, neither of which are solutions. (Figure El) Let us assume that in these 
regions, no solution exists, although our observations do not necessarily prove this. 
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Furthermore, let us adopt the following physical interpretation. By comparison, in 
the standard BCS theory of superconductivity, one has a finite temperature gap 
equation. As the temperature is raised, at the critical temperature there are no 
longer solutions to this gap equation; i.e. as far temperature is concerned, it is a 
bottom up approach. In contrast, the present formalism is a top down approach: 
as the temperature is lowered one reaches a critical temperature where solutions no 
longer exist. Let us interpret this as an instability toward formation of a new phase, 
or perhaps a cross-over to a different behaviour. In support of this interpretation, we 
mention the treatment of the unitary Bose gas within the present formalism|14l]. The 
gas undergoes a phase transition to a Bose-Einstein condensate at a critical value of 
fi/T; above this value there are no solutions to the pseudo-energy integral equation. 

It should be emphasized that the true nature of this 'phase' cannot be surmised 
from our thermodynamic approach alone; in addition one needs a bottom up ap- 
proach that contains information about the zero temperature ground state, such as a 
gap equation. Such a complementary bottom up approach is developed in 28[, where 
solutions to a BCS-like gap equation based on the attractive interactions described 
in section III are studied. The solutions are highly anisotropic, in that they vanish 
in the nodal directions, and are largest in the anti-nodal, and the critical T c ps 0.04 
found there for h = 0.15 is consistent with the critical temperatures estimated below. 

Figure M indicates the regions where solutions do not exist for positive hole doping 
h < 0.25. This figure is a contour plot of an interpolating function defined to be equal 
to 1 if there is a solution, and zero otherwise; the white region indicates the region of 
no solution, whereas in the light blue region, solutions exist. The boundary between 
the regions of existence and non-existence of solutions are the darkest curves, which 
are reasonably well delineated. Due to the 2-body approximation we have made in 
the thermodynamics, this boundary is not to be taken as precisely determined. The 
roughness of the boundaries we believe is a numerical artifact, mainly attributed to 
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not using a fine enough grid in the temperature and chemical potential variables. 
The dip around h = 0.15 we also believe to be an artifact since it disappears upon 
varying g and t' . In the range of doping 0.03 < h < 0.2 one sees a possible phase 
transition with critical temperatures ranging from < T c < .05. As explained in 
Section III, since this is the range of doping where the Fermi surface is intersecting 
the attractive band in the anti-nodal directions, we propose that this signifies an 
instability toward the formation of Cooper pairs, so that superconductivity may 
occur in the white regions. At hole doping h = 0.15, T c m 0.02. This is reasonable, 
since experimentally T Cimax /£ p» 0.025. 

The single quasi-particle energies correspond to e(k) + /i. In Figure [9] we plot this 
single particle energy as a function of temperature at optimal hole doping h = 0.15 
in the anti-nodal direction. One clearly sees a drop at T c>max . 

The higher T c 's up to 0.05 in the strongly underdoped region possibly signify the 
so-called pseudogap scale T*. There appears to be a small separation around h = .08, 
however this is less pronounced for other g, so it is not clear if this signifies anything. 
Some recent experimental results are very relevant to the issue 25|, |26[. Remarkably, 
it was found that the superconducting gap smoothly evolves into the pseudogap, 
i.e. they both seem to arise from the same underlying mechanism. In other words 
the gap is physically present even in regions of no superconductivity. This suggests 
that the T c 's in Figured] may all be arising from the same underlying phenomenon. 



This is consistent with the complimentary gap equation analysis in [28], where it 
was found that the gap extends and increases all the way to zero doping, as does the 
critical temperature scale in Figure |HJ Although not shown in Figure [8] at higher 
doping there is another region of no-solutions with a maximum T/t w 0.07. This 



could perhaps signify the temperature referred to as T co h in the 
crossover in the resistivity is observed from pocTtopcxT + T 2 



iterature, where a 



27]. 



We turn next to the thermodynamic functions, such as energy and entropy per 

21 



0.05 - 



0.04 



0.03 - 



0.02 - 



0.01 



0.00 - 




0.00 



0.05 



0.10 



0.15 



0.20 



0.25 



FIG. 8: Existence of solutions based on the iterative method. In the light regions there are 
no solutions to the integral equation for the pseudo-energy. The horizontal axis is the hole 
doping h, and the vertical axis is the temperature T. (See text for a detailed explanation.) 

particle. For hole densities in the vicinity of the boundaries shown in Figure [U our 
crude solution to the integral equation for the pseudo-energy is not smooth enough 
to reliably compute temperature derivatives numerically. However, at low density 
our numerical results are better behaved, and although of less interest physically 
for the cuprates, at least allow a comparison with previous literature. We therefore 
analyzed the thermodynamics in the overdoped region, with h = 0.8. The Fermi 
surface is shown in Figure [TO] In Figures [TTJ [12] and [131 we plot the energy and 
entropy per particle and specific heat as a function of temperature. Our results 
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FIG. 9: The single particle energy e(k) +/j, in the anti-nodal direction k = (0, ir) at optimal 
hole doping h = 0.15 as a function of temperature. 



for the entropy and specific heat are roughly consistent with the results in [8|, |9j], 
especially the results in [8j, which extend to low density; a detailed comparison is 
beyond reach since previous results are typically at higher temperatures. 
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FIG. 10: The Fermi surface for hole doping h = 0.8 at temperature T = 0.2. The white 
region has / = 1 and the darkest (purple) region / = 0. 

VI. CONCLUSIONS 



We presented an analytic treatment of the two-dimensional Hubbard model at 
finite chemical potential and temperature based on a new approach to statistical 
mechanics we recently developed ll|, |l2] , which is built upon the S-matrix. The ef- 
fective momentum-dependent coupling in this approach is the kernel G of an integral 
equation satisfied by the pseudo-energy e, which is built on the exact two-body S- 
matrix. We showed that there are regions of the Brillouin zone where the interactions 
are effectively attractive, for example, U/t > 6.4 for if ft = —0.3, even though the 
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FIG. 11: The energy per particle as a function of temperature for hole doping h = 0. 
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FIG. 12: Entropy per particle as a function of temperature for hole doping h = 0.8. 

bare model only has repulsive interactions, and this is essentially due to multi-loop 
quantum corrections. The next-to-nearest neighbor hopping coupling t' plays a sig- 
nificant role in determining this property, and our analysis suggests that for a fixed 
value of U/t, superconductivity may not exist for t' = 0. We emphasize that no ap- 
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FIG. 13: Specific heat per particle as a function of temperature for hole doping h = 0.8. 

proximations were made in obtaining these results; e.g. the existence of a threshold 
in U/t for the existence of certain attractive interactions, stands on its own and is 
independent of the subsequent approximations we made in the thermodynamics. 

We postulated that phase transitions occur where there are no solutions of the 
integral equation for the pseudo-energy. On the overdoped side, this phase sets in at 
hole doping h < 0.25. Our result for T c « 0.02 at h — 0.15 is in good agreement with 
experiments. We found that there is also evidence for transitions in the underdoped 
region, and we suggested this signifies the pseudogap. 

The main lesson of this work is that quantum loop corrections to scattering are 
perhaps the origin of the attractive interactions that lead to Cooper pairing near the 
Fermi surface. If this idea is correct, then in order to complete the picture one needs 
to derive a gap-equation that describes the structure of the ground state at zero 



temperature. Some preliminary attempts in this direction were taken in 28 



where 



solutions to a gap equation based on the attractive interactions described above are 
studied, and the critical temperatures found are consistent with the bottom down 
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approach described in this paper. Based on the detailed properties of the solutions 
to the gap equation studied there, i.e. its anisotropy and the existence of Fermi arcs 
in the nodal direction, it was suggested that the attractive mechanism in this paper 
may be responsible for the pseudogap rather than d-wave superconductivity. 
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